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Abstract 

Two  commonly  used  types  of  high-order-accuracy  element-based  schemes,  collocation- 
based  spectral  multidomain  penalty  methods  (SMPM)  and  nodal  discontinuous 
Galerkin  methods  (DGM),  are  compared  in  the  framework  of  the  inviscid  shal¬ 
low  water  equations.  Differences  and  similarities  in  formulation  are  identified, 
with  the  primary  difference  being  the  dissipative  term  in  the  Rusanov  form  of 
the  numerical  flux  for  the  DGM  that  provides  additional  numerical  stability; 
however,  it  should  be  emphasized  that  to  arrive  at  this  equivalence  between 
SMPM  and  DGM  requires  making  specific  choices  in  the  construction  of  both 
methods;  these  choices  are  addressed.  In  general,  both  methods  offer  a  mul¬ 
titude  of  choices  in  the  penalty  terms  used  to  introduce  boundary  conditions 
and  stabilize  the  numerical  solution.  The  resulting  specialized  class  of  SMPM 
and  DGM  are  then  applied  to  a  suite  of  six  commonly  considered  geophysical 
flow  test  cases,  three  linear  and  three  non-linear;  we  also  include  results  for 
a  classical  continuous  Galerkin  (i.e.,  spectral  element)  method  for  comparison. 

Both  the  analysis  and  numerical  experiments  show  that  the  SMPM  and  DGM 
are  essentially  identical;  both  methods  can  be  shown  to  be  equivalent  for  very 
special  choices  of  quadrature  rules  and  Riemann  solvers  in  the  DGM  along  with 
special  choices  in  the  type  of  penalty  term  in  the  SMPM.  Although  we  only 
focus  our  studies  on  the  inviscid  shallow  water  equations  the  results  presented 
should  be  applicable  to  other  systems  of  nonlinear  hyperbolic  equations  (such 
as  the  compressible  Euler  equations)  and  extendable  to  the  compressible  and 
incompressible  Navier-Stokes  equations,  where  viscous  terms  are  included. 
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1.  Introduction 


Geophysical  flows  exhibit  a  complex  structure  and  dynamics  over  a  broad 
range  of  scales  that  render  their  numerical  simulation  a  formidable  task  for 
state-of-the-art  computational  methods  and  resources.  Through  a  complex  in¬ 
terplay  between  the  earth’s  rotation,  ambient  stratification  and  the  constraining 
effects  of  lateral  and  vertical  boundaries,  flow  processes  in  geophysical  fluids 
commonly  exhibit  a  characteristic  horizontal  lengthscale  that  can  be  a  few  or¬ 
ders  of  magnitude  larger  than  its  vertical  counterpart  [1].  Hydrostatic  wave 
motions  occur  from  the  basin/planetary  scale  roughly  down  to  the  mesoscale. 
As  the  wave  scales  decreases,  non-linear  effects  become  significant  in  the  form 
of  internal/surface  bores  [2,  3].  At  wavelengths  of  O(lfcm),  the  waves  also  be¬ 
come  strongly  non- hydrostatic  [4].  The  fully  non-linear  and  non- hydrostatic 
waves  propagate  nearly  non-dissipatively  and  non-dispersively  over  long  dis¬ 
tances.  Turbulent  events,  driven  by  wave  breaking,  current-topography  inter¬ 
actions  and  other  mechanisms,  can  be  highly  localized  in  space  and  time  and 
span  a  broad  range  of  scales  within  their  region  of  occurrence.  Finally,  the  dis¬ 
sipative  effect  of  molecular  viscosity  is  only  felt  at  the  smallest,  O(lmm),  scales 
of  the  flow  field. 

As  a  result,  the  numerical  methods  used  in  the  investigation  of  geophysi¬ 
cal  flows  need  to  exhibit  a  number  of  preferred  features.  These  include:  a) 
front/wave  propagation  that  is  effectively  non-dissipative  and  non-dispersive, 
b)  minimum  artificial  dissipation  at  the  smallest  resolved  scales  to  enable  as 
broad  a  scale  separation  as  possible,  c)  efficient  resolution  of  localized  flow  fea¬ 
tures  and  complex  geometries  and  cl)  optimal  use  of  computational  resources. 
High-order  accurate  element-based  schemes  [5,  6]  are  particularly  appealing  in 
addressing  such  needs.  These  schemes  combine  the  exponential  convergence 
and  weak  artificial  dissipation  and  dispersion  of  standard  single-domain  spec¬ 
tral  methods  [7]  with  the  spatial  adaptivity  of  classical  finite  element /volume 
techniques  [8,  9].  Furthermore,  the  domain  decomposition  philosophy  inherent 
in  these  techniques  renders  them  highly  amenable  for  efficient  parallelization 
[!°]- 

On  account  of  the  inevitable  impossibility  of  capturing  the  full  range  of 
scales  intrinsic  to  a  highly  nonlinear,  and  steep,  front/wave  or  any  resulting 
localized  turbulent  event,  geophysical  flow  simulations  are  inherently  under¬ 
resolved.  Under-resolved  high-order  simulations  are  prone  towards,  often  catas¬ 
trophic,  numerical  instability  as  Gibbs  oscillations  are  compounded  by  aliasing 
driven  by  the  nonlinear  terms  in  the  governing  equations  [11].  In  high-order 
element-based  simulations,  these  numerical  instabilities  are  most  pronounced  at 
the  element  interfaces  when  strong  continuity  of  the  solution  is  enforced  across 
neighboring  elements  [12]  as  is  typically  done  in  continuous  Galerkin  methods. 

In  discontinuous  high-order  element-based  methods,  neighboring  subdomains 
carry  separate  values  of  the  solution  at  a  fixed  spatial  location  thereby  relax¬ 
ing  the  constraint  of  strong  continuity  of  the  solution  and  significantly  mit¬ 
igating  the  above  concerns  of  numerical  instability.  The  two  prevalent  cate¬ 
gories  of  such  methods  are  spectral  multidomain  methods  (with  and  without  a 
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penalty  scheme)  [13,  14,  15,  16,  17,  18,  19]  and  discontinuous  Galerkin  meth¬ 
ods  (DGM)  [20,  21,  22,  23,  6,  24,  25,  26].  Spectral  multidomain  methods, 
first  introduced  by  Orzag  [27]  for  elliptic  problems  were  originally  formulated 
with  a  strong  enforcement  of  continuity  of  the  solution  and  its  derivative  at  the 
subdomain  interfaces.  Subsequently,  Kopriva  [13]  extended  this  approach  to 
hyperbolic  problems,  where  the  interfacial  patching  was  implemented  with  an 
upwind  scheme  based  on  a  modified  method  of  characteristics.  This  approach 
was  further  refined  through  introducing  a  correction  method  based  combination 
of  characteristic  information  at  the  interfacial  points  [14].  In  the  framework  of 
the  compressible  Navier-Stokes  equations,  Kopriva  [15]  introduced  a  penalty  for¬ 
mulation  to  patch  subdomains  when  higher  (i.e. ,  2nd)  derivatives  were  present. 
The  formulation  used  in  our  work  follows  the  Spectral  Multidomain  Penalty 
Method  (SMPM)  presented  by  Hesthaven  [16]  and  expanded  upon  by  Don  [19] 
but  implemented,  to  our  knowledge  for  the  first  time,  to  the  shallow  water 
equations.  In  the  SMPM,  the  strong  interfacial  patching  conditions  are  re¬ 
placed  with  a  linear  combination  of  the  governing  equation  and  the  patching 
condition,  the  latter  multiplied  by  an  appropriately  chosen  penalty  coefficient. 
On  the  other  hand,  DGM  are  based  on  a  Galerkin  weighted  residual  formulation 
where  the  integration  is  performed  at  the  level  of  an  individual  element.  Since 
adjacent  elements  are  not  continuously  coupled,  as  is  the  case  with  finite  and 
spectral  elements,  interfacial  flux  integrals  do  not  vanish  and  are  represented  in 
the  form  of  an  appropriately  chosen  numerical  flux  that  preserves  consistency 
and  numerical  stability. 

SMPM  have  been  successfully  applied  to  high  Re  incompressible  stratified 
flow  process  studies  in  vertically  non-periodic  domains  such  as  internal  solitary 
wave-induced  bottom  boundary  layers,  turbulent  wakes  and  propagating  inter¬ 
nal  wave  packets  [28,  29,  30].  DGM  have  been  effectively  used  in  the  simulation 
of  the  shallow  water  equations  (SWE)  both  on  the  sphere  and  on  planar  but 
fully  unstructured  domains  [20,  21,  23,  25,  24]  and  for  compressible  atmospheric 
models  [22,  26]. 

However,  the  literature  exploring  the  similarities  and  differences  of  the  SMPM 
and  DGM  is  limited  to  the  recent  work  by  Gottlieb  and  Jung  [31]  who  consid¬ 
ered  the  modal  form  of  SMPM  and  DGM,  both  in  Galerkin  (integral)  formu¬ 
lation.  Focusing  on  one-dimensional  conservation  laws,  that  particular  study 
established  the  equivalence  between  the  two  techniques  for  a  specific  value  of 
the  penalty  coefficient  and  emphasized  the  additional  flexibility  of  the  penalty 
scheme  in  varying  the  value  of  this  coefficient  in  space  and  time  and  splitting  the 
advective  flux  at  the  subdomain  interfaces,  which  provided  for  greater  stability 
in  regions  of  strong  inhomogeneity  of  subdomain  thickness.  The  trade-offs  of 
accuracy  vs.  stability  as  a  function  of  the  penalty  coefficient  value  were  also 
examined  as  was  the  potential  of  the  coefficient  truncation  method  [32]  in  sup¬ 
pressing  rapid  error  growth  when  using  high-order  polynomials  in  the  penalty 
method.  Finally,  the  impact  of  inconsistent  evaluation  of  integrals  (exact  versus 
numerical  quadrature)  in  the  left  and  right-hand  sides  of  the  modal  Galerkin 
formulation  of  the  penalty  method  was  also  considered  in  the  framework  of  lin¬ 
ear  and  nonlinear  problems.  Note  that  both  the  coefficient  truncation  method 
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and  the  issues  with  integral  evaluation  are  restricted  to  the  modal  Galerkin  form 
of  the  SMPM. 

No  investigations  are  known  so  far  that  compare  the  collocation-based  SMPM 
and  the  nodal  Galerkin  formulation  of  the  DGM,  the  most  commonly  used  for¬ 
mulations  of  the  two  methods  which  this  paper  focuses  on.  Furthermore,  we  are 
unaware  of  any  comparison  of  the  two  methods  in  the  framework  of  a  system 
of  multi-dimensional  equations,  particularly  in  a  geophysical  context.  Such  a 
comparison  is  the  objective  of  the  present  paper.  The  platform  for  this  com¬ 
parison  are  the  SWEs  for  a  variety  of  reasons:  a)  the  relative  facility  of  their 
spatial  and  temporal  discretization  with  respect  to  more  complex  partial  dif¬ 
ferential  equations,  such  as  the  Navier-Stokes  equations,  b)  their  capability  for 
non-dissipative  propagation  of  highly  non-linear  waves,  which  renders  them  an 
ideal  experimentation  tool  for  testing  numerical  schemes  for  nonlinear  advec- 
tion,  the  primary  source  of  the  aliasing-driven  instabilities  mentioned  above  and 
c)  their  role  as  a  predictive  tool  of  ocean  wave  phenomena  for  the  purpose  of 
coastal  engineering  applications  [33]  and  tsunami  propagation  [34] .  We  specifi¬ 
cally  aim  to  compare  the  two  methods  in  terms  of  formulation  (with  a  focus  on 
subdomain  communication),  accuracy,  conservation  properties,  numerical  sta¬ 
bility  and  computational  cost  in  the  framework  of  specific  linear  and  non-linear 
test-cases. 

The  remainder  of  the  paper  is  organized  as  follows.  The  inviscid  SWE  are 
introduced  in  §2  along  with  their  representation  in  linear  and  quasi-linear  form. 
The  formulation  of  SMPM  and  DGM  is  presented  in  §3  along  with  an  overview 
of  the  accompanying  temporal  discretization.  SMPM  and  DGM  are  applied 
to  six  basic  test  cases  in  §4  followed  by  a  comparative  discussion  of  the  two 
methods  in  §5.  Conclusions  are  offered  in  §6. 

2.  Inviscid  Shallow  Water  Equations 

The  inviscid  shallow  water  equations  (SWE)  govern  the  behavior  of  a  fluid 
with  a  horizontal  extent  much  larger  than  its  depth,  and  are  derived  by  applying 
the  hydrostatic  approximation  to  the  incompressible  Navier-Stokes  equations 
[35] .  The  primitive  variable  formulation  of  the  SWE  is  given  by 
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where  tt,  v  are  the  horizontal  velocities,  H  is  the  mean  depth,  h  is  the  displace¬ 
ment  of  the  free  surface,  Z(u,  v )  is  the  external  forcing  and  g  is  the  gravitational 
constant. 
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2.1.  Conservative  form  of  the  SWE 

The  inviscid  shallow  water  equations  (equations  (1),(2)  and  (3))  can  also  be 
written  in  conservative  form: 


9q  *9F(q)  c>G(q) 

dt  dx  dy 

where  the  conservative  variables  q  are 


S(q), 


~  <f>~ 
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(4) 
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the  horizontal  and  vertical  fluxes  F(q)  and  G(q)  are  defined  as 
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and  the  source  terms  S(q)  are 
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In  Eqs.  (5),  (6)  and  (7),  <j>  =  gh  is  the  geopotential  height,  f  =  fo  +  P(y  ~  Um ) 
is  the  Coriolis  force,  tx  ,  tv  are  the  components  of  the  wind  stress,  p  is  the  fluid 
density,  and  7  is  a  bottom  friction  constant. 


2.2.  Linearized  SWE 

Assuming  a  mean  depth  much  larger  than  the  free  surface  elevation  (H  » 
h),  and  neglecting  the  nonlinear  terms  in  (4),  a  linearized  version  of  the  conser¬ 
vative  SWE  is  obtained.  The  modified  set  of  conservation  variables  is  defined 
as 
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where  $  =  gH  is  the  mean  depth  geopotential  height. 


2.3.  Quasilinear  form  of  the  SWE 

Using  the  chain  rule,  Eq.  (4)  can  be  rewritten  in  the  quasi-linear  form  [9,  36] 


dq  dF(q)  dq  dG(q)  dq  _ 

dt  <9q  dx  <9q  dy  tqJ 


where  A  and  B  are  the  flux  Jacobian  matrices,  that  can  be  decomposed  (via  an 
eigendecomposition  or  characteristic  decomposition)  as 

A  =  SaAaS^  (10) 

B  =  SbAbSjj1  (11) 

where  Aa  and  Ab  are  diagonal  matrices  containing  the  eigenvalues  of  A  and 
B,  and  Sa,  Sb  are  orthogonal  matrices  whose  columns  are  the  respective  eigen¬ 
vectors. 

The  positive  and  negative  flux  vectors  (F+,  F~,  G+,  G_)  are  defined  by 


F+  J 

[  SaA+S ~^dq 

(12) 

F  =J 

1  SaAaS  ~^dq 

(13) 

G+  = 

J 

1  SbA+S  ^dq 

(14) 

G~  = 

j 

1  SbAbS  ^dq 

(15) 

where  Aa  and  Ag  are  the  diagonal  matrices  composed  of  positive  and  negative 
eigenvalues  of  A  and  B,  respectively.  Based  on  the  above  decomposition,  the 
flux  vectors  have  the  properties 

Aa  =  A++A^  ->  F  =  F++F-  (16) 

Ab  =  A++Ab  ->  G  =  G+  +  G_.  (17) 

The  eigenvalue  matrices  and  flux  vectors  are  the  building  blocks  for  the  penalty 

formulation  of  the  SWE  via  SMPM,  and  for  the  definition  of  the  numerical  flux 

of  the  DGM  used  in  this  work  [20]. 

3.  Numerical  Methods 

3.1.  Spectral  Multidomain  Penalty  Method  (SMPM) 

The  SMPM  implemented  in  this  work  is  based  on  the  formulation  first  in¬ 
troduced  by  Hesthaven  [18]  and  further  refined  by  Don  et.al.  [19].  Specifically, 
this  SMPM  consists  of  a  multidomain  collocation  approach  based  on  discontin¬ 
uous  non-overlapping  rectangular  subdomains  that  are  connected  by  a  penalty 
term  that  ensures  stability  of  the  solution  by  imposing  weak  continuity  at  the 
subdomain  interfaces.  On  account  of  the  intrinsic  discontinuity  of  the  method 
and  the  critical  role  of  interfacial  patching,  the  formulation  of  the  SMPM  will  be 
presented  in  two  parts:  the  subdomain  interior  and  the  treatment  of  interfaces. 


6 


3.1.1.  Subdomain  Interior 

The  SMPM  is  based  on  a  collocation  approach  in  2D  quadrilateral  discon¬ 
tinuous  subdomains,  where,  within  each  subdomain,  any  function  q(x,  y,  t)  can 
be  approximated  by  using  N- th  order  Lagrange  interpolating  polynomials  on  a 
Gauss-Lobatto-Legendre  (GLL)  grid  [18]  as 

N  N 

q{x,y,t)  =  '^2'^2q(xi,yj,t)li(x)lj(y)  (18) 

£=0  j— 0 

where  q{xi,yj,t)  is  the  value  of  the  function  at  the  discrete  point  ( Xi,yj ),  and 
li(x),lj(y)  are  the  i— th  and  j— th  Lagrange  interpolating  polynomials  based  on 
the  GLL  nodes  in  the  x  and  y  directions,  respectively.  The  spatial  derivatives 
in  the  x— direction  in  the  global  coordinate  system  are  approximated  as 

dq(xi,yj,t)  dq(xi,yj,t)  d£ 

— Oi —  = — m — as  =  (19) 

s  fc= o 

where,  here,  we  assume  that  x  =  x(f)  and  f  =  f(x)  with  y  ^  y{x).  In  Eq.  (19) 
d£/dx,  represents  the  mapping  from  the  local  coordinate  system  £  €  [—1,1], 
given  by  the  GLL  points,  to  the  global  coordinate  system  x  £  R,  and  Dtj  is  the 
Legendre  spectral  differentiation  matrix,  that  is  computed  following  Costa  and 
Don  [37].  The  (/—derivative  is  approximated  in  a  similar  manner. 

3.1.2.  Interfacial  Treatment  and  Boundary  Conditions 

The  penalized  form  of  the  SWE  at  a  collocation  point  located  along  the 
boundaries  of  a  subdomain  requires  that  (see  reference  [19]  for  a  similar  for¬ 
mulation  of  the  compressible  Navier  Stokes  equations  for  chemically  reacting 
flow) 


=  S(q) 

+  r1Q(x)[F+(q)-  F+(q*)] 

+  t2  Q(x)  [F- (q)  -  F"  (q*)] 

+  t3  Q(x)  [G+  (q)  -  G+  (q*)] 

+  r4fi(x)[G-(q)-G-(q*)].  (20) 

In  (20),  Tj  (i  =  1,  •  •  •  ,4)  are  the  penalty  coefficients,  Q(x)  are  effectively  Dirac 
delta  functions  that  are  non-zero  only  at  the  interfaces  of  the  subdomain,  where 
the  penalty  terms  are  active,  and  F±(q),  G±(q),  F±(q*),  and  G±(q*)  represent 
the  positive  and  negative  fluxes  at  the  grid  points  on  the  particular  interfaces  of 
the  subdomain  (with  *  indicating  the  corresponding  point  on  the  neighboring 
interface)  on  the  subdomain  under  consideration.  In  a  general  sense,  the  penalty 
coefficients  can  be  viewed  as  weighting  factors  for  the  positive  and  negative 
fluxes  across  the  interfaces. 


<9q  dF(q)  <9G(q) 

dt  dx  dy 


7 


In  what  follows,  the  penalized  form  of  the  SWE  will  be  presented  for  the 
case  of  structured  quadrilateral  grids  with  rectangular  subdomains,  where  the 
treatment  for  vertical  interfaces  is  determined  by  the  horizontal  fluxes  dF/dx, 
and  for  the  horizontal  interfaces  by  the  vertical  fluxes  dG/dy.  Embedded  in  the 
penalty  coefficients  t,  (i  =  1,  •  •  •  ,4)  are  mapping  factors  to  enable  consistency 
in  units  between  the  different  terms  in  Eq.  (20). 


Vertical  interfaces.  Figure  1  presents  a  schematic  of  the  vertical  interface  be¬ 
tween  subdomains  /  and  II,  where  L  or  R  represent  any  collocation  point  at 
the  left  and  right  edges  of  the  interface. 


Figure  1:  Vertical  interface 


Based  on  (20),  the  penalized  form  of  the  SWE  for  a  point  located  at  the  left 
edge  of  the  interface  is 

dqL  dFL  dGL 


dt 


dx 


dy 


=  S(q  y 


+  nQL[(  F+)L-(F+)*] 

+  r2QL[(  F~)l-(F~)r].  (21) 


is 


Similarly,  for  a  point  along  the  right  edge  of  the  interface  the  penalized  form 


dqR  qFR  qqR 

dt  dx  dy 


S(q  f 

T5Qfl[(F+)fl-(F+)L] 

T6Qfl[(F-)*-(F-)L]. 


(22) 


In  Eq.  (22)  T5,r6  are  the  corresponding  penalty  coefficients  for  the  right  edge 
of  the  interface. 


Horizontal  interfaces.  Figure  2  presents  a  schematic  of  a  horizontal  interface 
between  subdomains  I  and  III.  In  this  case,  B  and  T  represent  the  collocation 
points  along  the  bottom  and  top  edges  of  the  interface.  The  penalized  equations 
for  a  point  located  at  the  bottom  edge  of  the  horizontal  interface  are 


dqB  dFB 
dt  +  dx 


dGB 

dy 


=  S(q)J 


+  r3QB[(G+)B-(G+)T] 
+  t4Qb[(  G")b-(G-)t] 


(23) 


B 

I 


Figure  2:  Horizontal  interface 


whereas  for  a  point  located  on  the 

<9qT  3Ft  dGT 
dt  dx  ch/ 


top  side  are 
=  S(q)T 


+  r7QT[(G+)T 
+  r8QT[(G^)T 


(G+)b] 

(G-)b], 


(24) 


In  Eq.  (24)  T7,t8  are  the  corresponding  penalty  coefficients  for  the  top  edge  of 
the  interface. 


The  approach  of  Don  et  al.  [19,  38]  for  a  one-dimensional  conservation  law 
can  be  readily  extended  to  the  penalized  equations  (21)- (24)  to  show  that  the 
penalty  scheme  formally  conserves  mass.  Moreover,  the  energy  of  the  system 
can  been  shown  to  be  bounded  by  its  initial  value  [19,  38]  if 


<  1,  2colt2  >  1 


2 lobt3  <  1, 
2lorT5  <  —  1, 
2 ujtT7  <  —  1, 
U)LTi  —  Ll1RT5  =  1, 
LOBT3  —  hJT  T 7  =  1, 


2coBT4  >  1 
2 lurt6  >  —  1 
2 ujtt8  >  —  1 

UJLT2  —  U)RTq  =  1 
LOBT4  —  U1TTs  =  1 


where  uL ,  u>B ,  uR  and  u>T  are  the  GLL  quadrature  weights  assigned  to  points 
along  the  left,  bottom,  right  and  top  interfaces,  respectively.  For  a  uniform 
order  of  polynomial  approximation,  N,  in  each  subdomain  a  single  value  of 
lu  =  2 /N(N  +  1)  can  be  used  instead. 
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Implementation  Issues.  In  this  work,  the  averaging  method  [19,  39]  is  imple¬ 
mented  such  that  the  penalty  coefficients  for  positive  and  negative  fluxes  (Eqs. 
(21)-(24))  at  the  sides  of  the  interfaces  are  taken  to  be  equal.  This  leads  to 


T~L  =  =  72 


tb  =t3  =  74 


tr-t5  =  76 


TT  —  T~7  =  Ts 


J_d£ 

2ix  dx 
1 

u)  Ax 
1  dr] 
2ix  dy 
1 

wA  y 
1  <9£ 

2ix  dx 

1 

ix  Ax 
1  dr] 
2 ix  dy 

1 

xAy 


(25) 


(26) 


(27) 


(28) 


where  d£/dx,dr)/dy  are  the  mapping  factors  for  the  penalty  terms  acting  on 
vertical  and  horizontal  interfaces  respectively  (see  Eqs.  (48)  and  (49)  ).  This 
approach  ensures  stability  of  the  penalty  scheme.  Moreover,  the  positive  and 
negative  fluxes  of  Eqs.  (16)  and  (17),  have  been  lumped  into  a  single  total  flux 
in  the  penalty  term. 

The  penalized  SWE  (  eqs.  (21)-  (24)  )  may  now  be  recast  accordingly  for 
each  possible  orientation  of  subdomain  interfaces: 


•  Vertical  interfaces 

—  Left  edge  of  the  interface 

dqL  dFL  dGL 

dt  dx  dy 

—  Right  edge  of  the  interface 

dqR  qFR  gGR 

dt  dx  dy 

•  Horizontal  interfaces 

—  Bottom  edge  of  the  interface 

dqB  dFB  dGB 
dt  dx  dy 


S(q)i+rLQL[FL-Ffl]  (29) 


S(q)R+TRQR[FR-FL }  (30) 


S(q)B  +tbQb[Gb  ~Gt]  (31) 
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Top  edge  of  the  interface 


dqT  dFT  dGT 
dt  dx  dy 


S(q)T  +  rTQT[GT-GB] 


(32) 


Note  that,  in  this  scheme,  unlike  Hesthaven  [18]  no  special  formulation  is  used 
for  the  corners,  which  are  simply  treated  as  points  that  belong  to  two  edges  of  the 
same  subdomain  orthogonal  to  each  other.  This  simplified  approach  is  found  to 
be  more  stable  than  the  theoretically  derived  one.  In  addition,  the  formulation 
of  the  penalty  term  is  the  same  form  used  by  Hesthaven  [17,  18],  Don  et  al. 
[19]  and  the  2nd  author  of  this  paper  [12].  Variations  of  this  formulation  are 
possible  and  a  particular  one,  involving  the  incorporation  of  dissipative  Rusanov 
flux-like  term,  is  examined  in  more  detail  in  §5.3. 


Compact  Representation  of  the  SMPM.  A  compact  form  of  representing  Eqs. 
(29)  -  (32)  is 

^  1=1 

where  n/e’^  is  the  outward  pointing  unit  vector  in  the  direction  from  control 
volume  e  to  l, 


with  As  =  (Ax,  Ay)  depending  on  the  orientation  of  the  subdomain  interfaces. 


3.2.  Discontinuous  Galerkin  Method  (DGM) 

The  discontinuous  Galerkin  (DG)  discretization  of  SWE  (4)  is  as  follows:  we 
begin  with  the  governing  equations  in  continuous  flux-form 

^+V-F(q)  =  S(q).  (34) 

Next  we  introduce  a  basis  function  expansion 

(iV+l)2 

Qn(x)  =  (35) 

i—1 

where  if  represents  the  basis  functions  of  order  N  and  qi  are  the  solution  vari¬ 
ables  at  specially  chosen  interpolation  points;  in  this  work  they  are  chosen  to 
be  the  Gauss-Legendre-Lobatto  (GLL)  points  in  order  to  make  the  comparison 
with  the  SMPM  more  relevant  and  because  we  have  used  these  points  in  pre¬ 
vious  DG  formulations  (e.g.,  [20,  22]).  Using  Eq.  (35)  we  can  now  construct 
approximations  for  the  remainder  of  the  spatial  terms  in  Eq.  (34).  For  example, 
we  can  now  represent  the  flux  tensor  as 

Fn  =  F(qN)  (36) 
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and  the  source  function  as 

SN  =  S(qN).  (37) 

Upon  defining  these  expansions,  we  can  then  substitute  them  into  Eq.  (34), 
multiply  the  equations  by  a  test  function,  and  integrate  to  obtain  the  element¬ 
wise  integral  problem:  find  qN  6  S(tte)Vip  €  <S(f2e)  on  each  element  f le  such 
that 

X  (cw  +  v'Fn)  dne  =  X  %kSN  dne  (38) 

where  S  is  the  finite-dimensional  space 

S  =  {ip  G  C2(fl)  :  1p\ne  G  P/v(^e)Vfte}  , 


PN  is  the  polynomial  space  of  order  N  defined  on  Qe  and  the  union  of  these  Ne 
elements  defines  the  global  domain,  i.e.,  f l  =  (J^i  ^e-  Next,  we  integrate  the 
divergence  term  by  parts  to  get 


da U) 

dne 

ot 


4 

£ 

i=i  • 


l/liTl 


(e,l)  _  p(e) 


N 


dTe 


[  V 

J  iie 


dfle 


dPLe 


(39) 


where  n^e’1^  is  the  outward  normal  vector  going  from  element  e  to  element  l 
that  defines  a  specific  edge  of  the  (in  this  specific  case)  quadrilateral  control 
volume.  Now,  since  the  solutions  are  discontinuous  across  element  boundaries 
then  it  becomes  critical  (in  order  to  construct  a  consistent  and  stable  numerical 
approximation  to  the  governing  continuous  equations)  to  choose  the  flux  tensor 
carefully.  To  resolve  this  inconsistency,  a  numerical  flux  is  introduced  that 
we  denote  by  The  simplest  choice  is  the  mean  value  between  the  two 

elements  claiming  the  same  interface 


p(*>0 

b  N 


+  F 


(O' 

N 


where  the  superscripts  e  and  l  represent  the  element  under  consideration  and 
the  side  (interface)  neighbor;  unfortunately  this  numerical  flux  is  not  the  best 
choice.  Another  easy  but  better  choice  is  the  local  Lax-Friedrichs  (or  Rusanov) 
flux  defined  as 


P(*,0  __  1 

*  N  ~  o 


p(e) 

N 


+  —  ddiss  | 


,(e,0 


(•In 


1 In' 


(e)^ 


(40) 


where  A  max  is  the  maximum  wave  speed  of  the  shallow  water  equations  (the 
maximum  eigenvalue  of  the  Jacobian  matrix  at  the  edge  l)  and  we  have  included 
the  switch  6(nss  that  controls  whether  the  dissipation  term  is  included.  With  a 
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specific  numerical  flux  defined,  the  DG  formulation  becomes 


da(e)) 

Vrt -  dne 

at 


4 

£ 

1=1  ■ 


tpiri 


(e,l)  .  p(*,0 


N 


dTe 


J  Q,e 


dfte 


ipiSff)  dfle 


(41) 


that  is  in  fact  the  weak  form  DGM.  Integrating  by  parts  one  more  time  yields 
the  following  mathematically  equivalent  system 


+  E  /  'il)in(e’l)  ■ 
1=1 


0  _  T?(e) 

*  N 


)  dTe+  [  dQe 

7  Jne 

=  [  ilnSff  dne  (42) 
Jfle 


which  is  the  strong  form  DGM  and  is  the  form  that  we  shall  use  to  compare  and 
contrast  with  the  SMPM  described  in  §3.1.  Next,  let  us  expand  the  terms  qN 
and  Sn  in  order  to  rewrite  Eq.  (42)  in  matrix- vector  form.  Expanding  these 
terms  in  Eq.  (42)  gives 


M 


(e) 


dq )' 


(e) 


dt 


£  «})r  (f£} 


=  M\fsf) 


(43) 


where  the  elemental  matrices  are  defined  as  follows: 


Ml 


[  i/jiijjj  dfle,  -D,)e- 

Jne 


[  dne,  M?) 

j  He 


(44) 


where  T  denotes  the  transpose  operator.  At  this  point  in  the  DG  formulation, 
we  have  to  introduce  numerical  quadrature  in  order  to  evaluate  the  integrals 
defined  in  Eq.(44)  in  the  following  way 


(Q+i)2 

Mif  =  £  ^k)\Jk)\Mxk)il’j(xk), 

k=  1 
(Q+1)2 

Dif  =  £  Uk)\Jk)\i’i(xk)'Vii>j{xk), 

k=  1 

(Q+1) 

Mij  =  £  ^k)\Jk)\^i(xk)i’j(xk)  (45) 

k= 1 

where  Q  is  the  number  of  quadrature  points  along  each  direction  of  the  quadri¬ 
lateral  element,  and  u  and  J  are  quadrature  weights  and  Jacobians,  respectively. 
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Using  GLL  points  for  both  interpolation  and  integration  we  obtain  the  fol¬ 
lowing  element  matrices 


Ml? 

(e)  I 

=  1 

■if 

DV 

(e)  I 

=  1 

jfWiiXi), 

M<‘> 

(/), 

=  “I  \. 

(46) 

where  <5  denotes  the  usual  Kronecker  delta  function.  Using  Eq.  (46)  in  Eq.  (43) 
and  dividing  by  the  mass  matrix  yields: 


dq 


(e) 


dt 


(V^ 


(xi)f  Ff 


=  S) 


(e) 


■X>,ro  cf'nf 

1=1 


(■ F f  -  ’°)  (47) 


where 


and 


Q 


(0 


1  if  i  is  on  the  edge  l 
0  otherwise 


T«  =  ^ 


P\jP\ 


(e) I  r(e)  I 


note  that  Eq.  (47)  is  quite  similar  to  Eq.  (33)  for  the  SMPM. 

Next,  we  need  to  simplify  the  penalty-like  term  that  we  have  called  r.  To 
do  so  requires  explicitly  stating  the  value  of  the  Jacobians  of  both  the  element 
and  edges.  For  the  sake  of  simplicity,  if  we  assume  that  £  =  £(x)  and  q  =  q (y), 
that  is,  that  the  computational  axes  are  aligned  exactly  with  the  physical  axes, 
then  we  can  write 


^  =  2(s  -  x0)  _  1 


V  = 


Ax 

2{y-yo) 

Ay 


-  1 


(48) 


where  xo,yo  is  the  left-bottom  most  point  on  each  element  and  Ax,  Ay  is  the 
length  of  the  element  along  the  x  and  y  directions,  respectively. 

This  mapping  yields  the  following  metric  terms 


_  J2_ 
dx  Ax 
dq  2 
dy  Ay 


(49) 


with  the  following  Jacobians 

1  r(e)  _  dx  dy  dx  dy  Ax  Ay 
dt;  dq  dq  dt;  4 
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and 


along  a  vertical  interface  (Left-Right  edge) 
along  a  horizontal  interface  (Top-Bottom  edge). 


\JW  1  = 


Ay 

2 

Ax 

2 


From  the  definition  of  these  metric  terms  we  can  see  that  the  penalty-like  term 
simplifies  to 


#  = 


-J27  along  a  vertical  interface  (Left-Right  edge) 

along  a  horizontal  interface  (Top-Bottom  edge) 


where  u>  =  ojq  =  lon  is  the  value  of  the  quadrature  weight  at  the  beginning  or 
end  point  (they  are  equal  by  symmetry).  Introducing  the  DGM  numerical  flux 
given  in  Eq.  (40)  into  Eq.  (47)  yields 


dq 


(e) 


dt 


+  (V$j(xi))T  F<f> 


=  S[e) 


EEs: 


p)n(i)n(e,0 


r>(e) 


-  F®  -  <5diss|A. 


diss  | ^max 


1=1 


„■>( 


(0  (e) 

Q  i  ~  <? 


(50) 


where 

^  _  T  1 

2  ccAs 

and  As  =  (Ax,  Ay)  depending  in  which  direction  the  interface  is  oriented.  At 
this  point,  we  have  not  made  too  many  sacrifices  or  simplifications  in  deriving 
Eq.  (50).  This  equation  is  in  fact  a  valid  DGM  representation  of  the  shallow 
water  equations  with  only  the  very  slight  assumptions  that: 

1.  The  computational  coordinates  (£,77)  are  aligned  with  the  physical  coor¬ 
dinates  ( x,y ). 

2.  Co-located  interpolation  and  integration  points  are  used.  The  fact  that 
we  have  chosen  these  points  to  be  the  GLL  points  results  in  inexact  inte¬ 
gration. 

3.  The  numerical  flux  used  is  the  simple  Rusanov  flux. 

Taking  the  special  case  Sdiss  =  0,  that  is,  no  dissipation  in  the  flux  term,  yields 


j  (e) 

dq) 

dt 


+  (vyj(xi))T  Ff 


s[e) 


+  E^W 

1=1 


(51) 


which  is  identical  to  the  SMPM  representation  given  in  Eq.  (33).  Eq.  (51) 
shows  that  another  way  of  viewing  the  penalty  term  is  as  an  extra  differencing 
term  (as  is  evident  by  the  -g-  term  in  r  and  A F  in  the  numerator)  that  considers 
the  information  from  the  neighboring  elements,  which  is  in  fact  what  we  mean 
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s 

&ik 

fiik 

i 

1/2 

4 

0 

1 

0 

1/2 

2/3 

0 

1/3 

0 

0 

1/6 

0 

0 

0 

1 

0 

0 

0 

1/2 

Table  I:  Coefficients  for  the  third  order  -  four  stage  SSP-RK(34)  method 


by  the  usual  term  flux.  In  §4  we  use  Eq.  (50)  with  and  without  the  dissipation 
term  to  compare  the  SMPM  with  the  DGM.  We  now  turn  our  discussion  to  the 
time-integrator  we  use  to  advance  the  SMPM  and  DGM  solutions  forward  in 
time. 

3.3.  Temporal  Discretization 

To  retain  the  high-order  accuracy  of  the  SMPM  and  the  DGM,  a  high-order 
time  advancement  scheme  is  needed.  The  explicit  strongly  stability  preserving 
Runge-Kutta  (SSP-RK)  method  [40,  41]  is  implemented  for  both  approaches. 
Consider  the  following  initial  value  problem 

I  =  R^-  <52> 

The  prediction  at  the  time  n  +  1  is  based  on  the  existing  solution  at  the  time  n 
and  the  forcing  terms  R(q).  The  scheme  can  be  written  as  [41] 


9<°) 

=  qn 

(53) 

q(i) 

i—  1 

=  5l(Wfc)+  At/3ikR{qWj)  , 

k— 0 

*  =  1,2,.. 

-,s  (54) 

,(«+!) 

=  g(s) 

(55) 

where  s  are  the  number  of  stages  of  the  SSP-RK  approach,  and  ffk  are 
constant  coefficients  given  in  Table  I  [41],  and  At  is  the  size  of  the  time  step  at 
a  specific  time. 

4.  Test  cases:  Description  and  Results 

Six  test  cases  are  examined  to  compare  the  performance  of  the  SMPM  and 
DGM  in  terms  of  accuracy,  dynamic  stability,  robustness  and  conservation  prop¬ 
erties:  three  linear  (standing  wave,  Kelvin  wave,  and  Stonnnel  problem),  where 
accuracy  can  be  evaluated  through  the  availability  of  analytic  solutions,  and 
three  non-linear  (nonlinear  Stommel,  equatorial  Rossby  wave,  and  Riemann 
problem)  that  provide  a  platform  for  assessing  the  dynamic  stability  and  ro¬ 
bustness  of  the  methods.  In  addition,  results  obtained  with  the  spectral  element 
method  (SEM)  [42,  43,  44]  are  included  to  compare,  for  each  case,  the  behavior 
of  a  continuous  method  with  a  discontinuous  element-based  approach.  For  the 
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linear  cases  an  additional  error  analysis  based  on  the  normalized  L ^  and  L2 
norms  of  the  error  is  performed. 


The  normalized  L ^  and  L2  error  norms  are  defined  as 


£=  = 


I \Hl2  = 


maxxe n{hexact  -  h) 
maxxg  q /i.exac£ 

I  ffoi^exact  h)2dfl 
Jq  h exact 


(56) 

(57) 


The  mass  (M)  and  energy  (E)  of  the  system  are  measured  in  the  following  way 

M  =  I  4>dn  (58) 

Jn 

E  =  [  [c/){u2  +  v 2)  +  (f> 2]  dn.  (59) 

Jn 


The  metric  for  assessing  mass  and  energy  conservation  is  the  respective  relative 
error,  defined  with  respect  to  the  corresponding  initial  values  of  M  and  E.  It 
is  computed  as 


R\i 


Mt  -  M0 
M0 


Re 


Et  —  Ep 

Eo 


(60) 


where  Rm  and  Re  are  the  relative  errors  in  mass  and  energy,  and  Mo,  Eo,  Mt,  Et 
are  the  corresponding  values  for  mass  and  energy  at  the  initial  and  final  times 
of  the  simulation,  respectively.  For  each  test  case,  it  is  specified  explicitly  if 
mass  and  energy  are  lost  or  generated  by  the  end  of  simulation. 

For  all  simulations  no  boundary  conditions  are  applied  to  the  continuity 
equation.  For  the  momentum  equation  no-flux  (i.e. ,  reflecting)  boundary  con¬ 
ditions  are  applied  along  all  four  walls  of  the  basins;  for  the  SEM  and  SMPM 
methods  this  is  accomplished  via  strong  homogeneous  Dirichlet  boundary  con¬ 
ditions  whereas  for  the  DGM  they  are  satisfied  in  a  weak  sense. 

To  compute  the  Courant  number  a  high-order  cell  technique  is  used,  where 
the  cells  are  defined  based  on  the  GLL  points  on  each  subdomain.  A  mean 
velocity  and  geopotential  height  is  defined  at  the  center  of  each  cell  [23] .  With 
these  considerations,  the  Courant  number  is  defined  as 

Courant  Number  =  max 

where  At  is  the  size  of  the  time  step,  U  is  the  mean  velocity  magnitude  at  the 
cell,  <j>  is  the  average  geopotential  height  in  the  cell  and  As  =  \J Ax2  +  Ay 2 
is  the  grid  spacing.  For  SMPM  and  DGM,  the  maximum  Courant  number  At 
that  ensures  stability  of  the  numerical  simulations  0.5  (Courant  Number  <  0.5). 
The  equivalent  value  for  SEM  is  1.  As  specified  in  the  relevant  sections,  two 


^ At(U  + 
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test  cases  (standing  and  Kelvin  wave)  are  run  with  a  significantly  smaller  time 
step  to  prevent  the  time-stepping  error  from  dominating  the  error  associated 
with  the  spatial  discretization.  Nonetheless,  as  the  conservation  properties  of 
the  SMPM  are  negatively  impacted  by  a  linearly  growing  loss  of  mass  which  is 
of  order  machine  epsilon  at  each  time  step,  all  other  test  cases  are  run  with  a 
time  step  that  is  80%  the  maximum  time  step  associated  with  Courant  number 
limits  indicated  above.  We  refer  the  reader  to  §5.1  and  5.2  for  further  discus¬ 
sion  on  time-stepping  error  and  the  impact  of  time  step  on  the  conservation 
properties  of  the  spatial  discretization  methods  under  consideration.  The  de¬ 
gree  of  polynomial  approximation  is  varied  from  N  =  4  to  20.  The  number  of 
subdomains  is  also  varied  within  a  range  dependent  on  the  geometry  of  each 
case,  and  the  SSP-RK34  method  defined  previously  is  used  to  advance  in  time 
the  simulations. 


4-1.  Linear  Problems 

In  this  section,  we  compare  the  three  methods  quantitatively  using  linear 
test  cases  that  have  analytic  solutions. 


4-1.1.  Linear  Standing  Wave 

This  case  represents  the  evolution  in  time  of  a  wave  driven  only  by  gravita¬ 
tional  effects  (S  =  0)  through  an  initial  perturbation  of  the  free  surface.  From 
references  [23,  45],  the  analytic  solution  for  this  case  is  given  by 


h(x,  y,  t) 
u(x,y,t) 

v(x,y,t) 


cos  (7 xx)  cos  (7 xy)  cos  (7 Ra/2) 
1 

75 

1 

75 


sin  (ttx)  sin  (7 xy)  sin  (ntV2) 
cos  (nx)  sin  (7 xy)  sin  (nt\/2) 


(61) 

(62) 

(63) 


with  (, x,y )  €  [0, 1]  x  [0, 1]. 

The  simulations  are  run  for  t  £  [0,  0.5].  Figure  3  shows  results  for  SMPM, 
DGM  and  SEM  simulations  for  a  fixed  number  of  subdomains  and  variable 
order  of  polynomial  approximation  N.  A  time  step  which  is  l/50th  of  that 
associated  with  a  Courant  number  value  of  0.4  is  used,  to  make  time-stepping 
errors  sufficiently  small.  The  results  are  indistinguishable  if  an  even  smaller  time 
step  is  employed.  Exponential  convergence  of  the  error  norms  for  free  surface 
elevation  and  horizontal  velocity  is  attained  for  each  method  for  polynomial 
degree  less  or  equal  than  N  =  8.  At  higher  values  of  N,  the  convergence  rate  is 
finally  reduced,  reaching  a  plateau  of  the  order  of  O(10-12),  the  cause  of  which 
we  are  unable  to  determine.  The  Galerkin  based  methods  (i.e.  DGM,  SEM) 
conserve  mass  up  to  machine  precision.  The  SMPM  mass  cumulatively  loses 
mass  over  time.  All  three  methods  show  improved  energy  conservation  with 
increasing  N  with  the  relative  error  reaching  a  value  of  0(  10-12)  at  IV  =  8.  An 
interpretation  for  the  performance  of  the  SMPM  in  terms  of  mass  conservation 
is  offered  in  §5.1. 
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Figure  3:  Analysis  of  the  standing  wave  (5x5  subdomains)  at  t  =  0.5  seconds  for  a  varying 
number  of  GLL  points,  a)  L2  normalized  relative  error  in  the  free  surface  elevation  h.  b)  L2 
normalized  relative  error  in  u  velocity,  c)  Relative  error  in  mass,  d)  Relative  error  in  energy. 

4-1-2.  Linear  Kelvin  Wave 

The  equatorial  Kelvin  wave  is  a  low  amplitude  non-dispersive  wave  trapped 
in  the  vicinity  of  the  equator.  It  is  driven  by  rotational  and  gravitational  effects 
through  an  initial  perturbation  of  the  free  surface.  The  analytic  solution  for 
this  case  [23,  46]  is 

h(x,y,t )  =  1  +  exp  ^  exp  ^  +  ^ ^  (64) 

u(x,y,t)  =  exP  (~y)  exP  ^  +  2 — — )  (65) 

v(x,y,t)  =  0  (66) 

for  f0  =  0,  P  =  1  and  (x,  y)  €  [-20,  20]  x  [-10, 10]. 

Simulations  are  run  for  t  £  [0,5].  Figure  4  shows  results  for  this  case  for 
a  domain  discretized  with  20  x  10  elements  and  a  varying  value  of  N.  As 
with  the  Standing  wave,  here  the  time  step  is  l/50th  that  associated  with  a 
Courant  number  value  of  0.4.  No  further  reduction  in  time  step  was  required  to 
make  time-stepping  errors  sufficiently  small.  The  behavior  of  the  error  norms  is 
similar  to  that  observed  for  the  linear  standing  wave:  exponential  convergence 
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is  observed  for  all  the  three  methods.  DGM  and  SEM  conserve  mass  up  to 
machine  precision.  On  the  contrary,  SMPM  again  shows  a  loss  of  mass,  which, 
in  the  end  of  simulations,  is  up  to  one  order  of  magnitude  larger  than  the 
value  computed  for  DGM  and  SEM.  The  trend  in  relative  error  of  total  energy 
conserved  is  comparable  to  that  observed  for  the  linear  standing  wave  in  Fig. 
3.  Improved  energy  conservation  occurs  with  increasing  N  with  a  relative  error 
value  of  O(10-13)  observed  for  N  =  20. 


b.  u-velocity  normalized  L  error  norm 


Figure  4:  Kelvin  wave  results  for  20  X  10  subdomains  at  t  =  5.  Panels  (a)  through  (d)  show 
the  same  quantities  with  Fig.  3. 


4-1.3.  Linear  Stommel  Problem 

This  problem  [47]  also  known  as  westward  intensification  of  wind-driven 
ocean  currents,  represents  the  steady  balance  between  rotation,  gravity,  friction 
and  wind  stress  in  a  square  ocean  basin.  A  sinusoidal  wind  stress  forces  an  un¬ 
perturbed  free  surface  generating  a  small  amplitude  wave  moving  westward  due 
to  the  Coriolis  force  that  is  compensated  by  bottom  friction  and  gravitational 
effects  and,  eventually,  reaches  steady  state.  The  analytic  solution  used  for  this 
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case  is  [23] 


where 


h(x,y,t) 


u(x,y,t) 
v(x,  y,  t) 


XlX  +  0-eX2X)-^-  cos 
A2  l 


\2  fny\ 

)  cos  It) 

\  7T  , 

+f  sin  (CieAl*  +  C2ex*x  +  C3) 
—  (CieXlX  +  C2ex*x  +  C3) y  cos 
(CiAieAl"  +  C2A2eA2a:)sin(^) 


Ci 

c2 

c3 


c3 


-Cl 

rl 

7T7 


1  -  eX21 

eX2l  _ gAi / 

1  -  eAli 


D  A2  l  pXll 


(67) 

(68) 

(69) 

(70) 

(71) 

(72) 


For  the  case  presented  here,  /o  =  lx  10~4,  /3  =  1  x  10-11,  7  =  1  x  10-6,  g  =  10, 
p  =  1000,  r  =  0.2,  F/o  =  1000,  and  (. x,y )  €  [0, 1  x  106]  x  [0, 1  x  106].  Note  that 
the  solution  is  symmetric  with  respect  to  the  y  axis. 


a.  SMPM  b.  DGM  c.  SEM 


Figure  5:  Free  surface  elevation  computed  by  all  three  methods  for  the  linear  Stommel  problem 
for  5x5  subdomains  and  =  12  at  t  =  400  days 


Simulations  are  run  until  the  solution  is  close  to  the  steady  state  (i.e.  t  =  320 
days),  and  the  structure  of  the  steady  state  flow  field,  displaying  the  expected 
symmetry  around  the  horizontal  axis  at  z  =  5  x  105,  is  shown  in  Fig.  5  for 
all  three  methods.  Figure  6  shows  the  error  norm  convergence  curves  for  the 
case  of  a  5  x  5  mesh  for  solutions  obtained  with  different  values  of  N.  For  all 
three  methods,  the  error  in  the  free  surface  displacement  shows  an  exponential 
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Figure  6:  Linear  Stommel  problem  results  for  5x5  subdomains  at  t  =  320  days,  a)  Con¬ 
vergence  plot  for  the  L2  normalized  relative  error  in  the  free  surface  elevation  h.  b)  Relative 
error  in  mass  conservation. 


convergence  similar  to  the  previous  two  linear  cases  for  up  to  N  =  8,  beyond 
which  the  error  norms  level  off  to  a  constant  value.  This  plateau  is  reached 
because  an  exact  steady  is  almost  never  attained  in  practice,  as  simulations 
are  dominated  by  slowly-decaying,  weak-amplitude  basin-scale  modes,  with  the 
decay  time  of  the  gravest,  longest- wavelength,  mode  reaching  60  years  [48]. 
Mass  is  conserved  up  to  machine  precision  by  DGM  and  SEM,  whereas  SMPM 
shows  a  loss  of  total  mass  up  to  three  orders  of  magnitude  larger  than  DGM 
and  SEM. 

4-2.  Nonlinear  Problems 

In  this  section,  we  compare  the  three  methods  qualitatively  using  nonlinear 
test  cases  that,  unfortunately,  do  not  have  analytic  solutions.  Instead,  we  use 
the  conservation  of  mass  and  energy  to  compare  the  methods.  All  three  models 
formally  should  conserve  mass  but  are  not  guaranteed  to  conserve  energy.  It  is 
possible  to  conserve  energy  (at  least  up  to  the  time-truncation  error)  but  this 
requires  slight  modifications  to  the  discrete  operators  that  we  will  not  pursue 
in  this  work. 

4.2.1.  Nonlinear  Rossby  Soliton 

This  case  considers  an  equatorial  non-linear  Rossby  wave  of  weak  amplitude, 
driven  by  gravity  and  rotational  forces.  It  is  initialized  by  a  Gaussian-like  per¬ 
turbation  in  the  free  surface  elevation.  An  approximate  asymptotic  solution 
of  the  system  of  Korteweg-DeVries  equations  resulting  from  the  SWE  through 
application  of  the  method  of  multiple  scales  is  obtained  for  this  problem  in 
[49].  Although  this  first  order  solution  does  not  provide  a  reference  to  assess 
the  convergence  rate  of  the  numerically  computed  solution  for  the  SEM,  DGM, 
and  SMPM,  it  is  used  to  compare  associated  phase  speed  and  solution  structure 
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with  the  corresponding  estimates  computed  by  the  three  numerical  methods. 
For  this  case  (x,  y)  £  [—24,  24]  x  [—8, 8],  g  =  1,  and  the  Coriolis  force  f(y)  =  y. 

Simulations  are  run  for  t  £  [0,40].  All  three  methods  accurately  reproduce 
the  free  surface/ velocity  structure  of  the  soliton  and  its  propagation  at  a  con¬ 
stant  phase  speed  equal  to  the  analytically  predicted  value.  The  structure  of 
the  free  surface  elevation  field  at  the  end  of  the  simulation,  with  its  character¬ 
istic  two-lobe  structure,  as  computed  by  all  three  methods  is  shown  in  Fig.  7. 
Figure  8  shows  results  for  mass  and  energy  conservation  for  24  x  8  subdomains, 
and  varying  N,  which  are  similar  to  their  counterparts  obtained  for  the  linear 
cases.  The  SMPM  is  subject  to  a  decrease  in  mass  when  the  polynomial  order 
increases.  The  DGM  conserves  mass  up  to  machine  precision,  with  the  SEM 
offering  comparable  performance.  The  SEM  and  SMPM  are  the  most  and  least 
energy  conserving,  respectively.  As  discussed  in  §5.1,  the  energy  conservation 
properties  of  the  DGM  are  highly  dependent  on  the  formulation  of  the  numerical 
flux  and  the  use  of  spectral  filtering  (see  Fig.  12). 


Figure  7:  Qualitative  comparison  of  the  Non-linear  Rossby  wave  results  with  24  x  8  subdo¬ 
mains,  N  =  12,  and  at  time  t  =  40. 


4-2.2.  Nonlinear  Stommel  Problem 

The  same  configuration  (forcing  parameters,  dimensions  of  the  physical  do¬ 
main,  and  boundary  conditions)  is  used  as  in  the  linear  Stommel  problem.  How¬ 
ever,  the  fully  nonlinear  set  of  Eqs.  (4)  are  now  solved.  In  this  case,  a  shift  of 
the  gyre  toward  the  northwest  part  of  the  basin  is  expected  due  to  the  effect  of 
the  nonlinear  terms. 

Figure  9  shows  the  steady  state  results,  for  a  domain  with  5x5  subdomains. 
Similar  trends  are  observed  for  all  three  methods.  Note  that  in  this  particu¬ 
lar  case,  the  differences  in  subdomain  interface  treatment  between  SMPM  and 
DGM  give  rise  to  challenges  of  numerical  stability  for  the  former,  when  values 
of  polynomial  degree  N  are  used.  In  the  SMPM,  when  5x5  subdomains  are 
used  and  N  >  12,  weak  spurious  oscillations  develop  in  the  top  left  corner  of 
the  domain  and  intensify,  as  time  advances,  eventually  forcing  a  catastrophic 
blow-up  of  the  solution.  As  a  counter-measure,  a  16-th  order  Boyd-Vandeven 
filter  [50]  is  used,  which  attenuates  only  the  very  highest  modes  of  the  solution, 
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Figure  8:  Non-linear  Rossby  wave  results  for  24  X  8  subdomains  at  f  =  40.  a)  Relative  error 
in  mass,  b)  Relative  error  in  energy. 

a.  SMPM  b.  DGM  c.  SEM 


Figure  9:  Nonlinear  Stommel  problem  for  5x5  subdomains  and  .V  =  1 2  at  t  =  400  days 


to  suppress  these  oscillations.  This  problem  does  not  occur  for  the  DGM,  as  the 
spurious  oscillations  are  damped  by  the  dissipative  term  Sdiss  =  1  hr  the  numer¬ 
ical  flux.  The  sensitivity  of  the  DGM  and  SMPM  to  the  presence  of  dissipative 
terms  is  examined  in  greater  detail  in  §5.3.  Figure  10  shows  the  behavior  of  the 
relative  error  in  mass  as  a  function  of  N,  which  is  similar  to  what  is  observed 
for  the  corresponding  linear  problem  (Fig.  6)  .  Results  are  restricted  to  N  <  8, 
as  high-order  polynomial  approximations  require  the  use  of  a  spectral  filter  to 
preserve  stability. 


4-2.3.  Nonlinear  Riemann  Problem 

This  modification  of  the  circular  dam  break  problem  [51]  is  considered  as  a 
platform  to  assess  the  performance  of  the  three  methods  in  simulating  strongly 
nonlinear  flows,  i.e.  flow  fields  with  distinct  sharp  spatial  gradients.  The  initial 
condition,  a  Gaussian  bump  (used  instead  of  a  cylindrical  step  function),  is 
characterized  by  such  a  sharp  gradient  and  has  free  surface  and  velocity  fields 


24 


10“ 

10' 

ffi  10" 

E 

§  10" 


10" 

10" 


Figure  10:  Relative  error  in  mass  as  a  function  of  polynomial  order  for  the  Nonlinear  Stommel 
problem.  5x5  subdomains  at  t  =  360  days. 

given  by: 

k(x,vM  -  H+Aex p(-(— o^-^)  (73, 

u(x,y,t0)  =  0 

v{x,y,t0)  =  0 

(74) 

where  ( x,y )  G  [0,1]  x  [0,1],  g  =  9.8,  H  =  1,  A  =  0.2,  Xq  =  yo  =  0.5,  and 
a  =  0.05.  The  flow  is  driven  by  gravity  as  in  the  standing  wave  problem.  Sim¬ 
ulations  are  run  for  t  £  [0,0.2],  i.e. ,  up  to  a  short  time  after  the  first  reflection 
of  the  initial  wave  from  the  domain  boundaries  where  reflecting  boundary  con¬ 
ditions  are  applied. 

Figure  11  shows  results  for  conservation  properties  in  the  case  of  a  5  x  5  sub- 
domains.  In  terms  of  mass  conservation,  it  is  difficult  to  discern  which  method 
offers  superior  performance.  The  energy  conservation  properties  of  each  method 
improve  with  increasing  N.  At  a  given  value  of  N ,  the  DGM  is  found  to  produce 
a  slightly  larger  relative  error  in  terms  of  the  total  final  energy.  Note  that  for 
the  time  for  which  the  simulations  were  run,  no  filtering  was  needed  to  preserve 
numerical  stability  at  all  values  of  N  and  subdomain  thicknesses  considered. 
Nevertheless,  the  smoothness  of  the  solution  is  damaged  at  later  times,  as  weak 
spurious  wiggles  emerge.  As  in  the  case  of  the  non-linear  Stommel  problem,  in 
the  DGM,  the  dissipation  term  in  the  Rusanov  flux  stabilizes  the  solution  while 
keeping  it  free  of  spurious  oscillation,  although  somewhat  adversely  impacting 
the  energy  conservation  properties  of  the  method.  The  role  of  spectral  filtering 
and  dissipative  terms  on  the  conservation  properties  for  the  DGM  is  further 
discussed  in  §5.1. 


Mass  conservation  (m2/s2) 


2  3  4  5  6  7  8 

Polynomial  degree  (N) 
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Figure  11:  Nonlinear  Riemann  problem  for  5x5  subdomains  at  t  =  0.2.  Panels  (a)  and  (b) 
are  the  same  as  Fig.  8 


5.  Discussion 

5.1.  Mass  and  Energy  Conservation 

All  three  methods  are  found  to  have  very  good  conservation  properties,  a 
direct  result  of  their  formulation,  see  e.g.  [19]  for  SMPM,  [22]  for  DGM,  and 
[52]  for  SEM.  The  DGM  conserves  mass  up  to  machine  precision.  The  SMPM 
is  found  to  lose  mass  over  long  model  times  with  the  corresponding  relative 
error  as  much  as  four  orders  of  magnitude  larger  than  that  for  the  DGM.  This 
error  increases  with  number  of  time  steps.  Such  observations  might  initially 
seem  perplexing,  given  the  analytical  demonstration  of  Don  et  al.  [19]  that 
the  averaging  method-based  penalty  scheme  is  conservative.  For  all  SMPM- 
driven  test-cases  we  have  found  that  the  mass  loss  (not  shown  here)  is  a  linear 
function  of  time,  with  a  decay  rate  that  is  of  the  order  of  machine  epsilon.  The 
linear  Stommel  problem  has  a  total  mass  loss  that  reaches  values  of  10“ 10  at 
higher  N,  a  value  even  higher  than  that  observed  for  the  standing  and  Kelvin 
wave  test-cases  where  l/50th  the  maximum  time  step  is  used.  This  difference 
is  simply  because  106  time  steps  are  required  for  the  linear  Stommel  problem 
to  reach  steady-state.  Consequently,  we  attribute  the  observed  loss  of  mass  to 
an  accumulation  of  round-off  error. 

The  energy  conservation  properties  of  all  three  methods  improve  with  in¬ 
creasing  N,  although  both  SMPM  and  DGM  are  found  to  be  inferior  in  this 
regard  to  the  SEM.  Note  that  in  simulations  where  no  energy  sink  terms  (such 
as  bottom  friction  in  the  Stommel  problems)  are  present,  the  performance  of 
the  discontinuous  techniques  in  terms  of  energy  conservation  can  be  strongly 
influenced  by  spectral  filtering  and  the  structure  of  the  numerical  flux  terms, 
such  as  the  dissipative  term  used  within  the  Rusanov  flux.  Figure  12  shows 
the  differences  in  conservation  of  mass  and  energy  in  the  DGM,  for  the  Rie¬ 
mann  problem,  when  spectral  filtering,  through  a  lOth-order  Boyd-Vandeven 
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filter  [50],  is  added  to  the  simulation  or  the  dissipation  term  is  neglected  in 
the  numerical  flux.  The  absence  of  both  the  dissipative  term  in  the  numerical 
flux  and  spectral  filtering  provides  for  the  best  energy  conservation  properties, 
although  such  behavior  does  not  necessarily  guarantee  a  smooth  and  stable  so¬ 
lution  for  such  a  strongly  nonlinear  problem. 


Figure  12:  Comparison  of  conservation  properties  of  the  DGM  for  the  Riemann  problem. 
Results  for  5  x  5  at  t  =  0.2.  (a)  Mass  conservation,  (b)  Energy  conservation. 


5.2.  Effect  of  time  step  on  convergence  and  conservation  properties 

For  the  purpose  of  demonstrating  that  the  temporal  discretization  error  does 
not  dominate  over  the  spatial  error,  we  now  perform  an  analysis  of  the  effect  of 
time  step,  At,  size,  on  the  convergence  and  conservation  properties  of  each  of 
the  three  methods.  The  base  time  step  corresponds  to  that  associated  with  a 
simulation  with  Courant  Number  of  0.4.  At  is  then  progressively  decreased  by 
a  factor  of  2, 10  and  50  (denoted  by  D 2,  D 10,  D50  respectively).  In  figure  13  the 
convergence  plots  for  the  free  surface  elevation  h  of  the  Standing  wave  test  case 
are  presented  for  all  three  methods.  For  a  given  TV,  the  increase  in  accuracy  of 
all  three  methods  is  visible  as  At  is  decreased.  Once  a  factor  of  50  reduction  is 
reached  exponential  convergence  is  obtained  until  N  =  8. 

The  same  exercise  has  been  performed  to  assess  the  role  of  time-step  on  mass 
and  energy  conservation  in  all  three  methods.  The  results  show  (see  figures  14 
and  15)  that  the  SMPM  mass  loss  increases  with  decreasing  At.  This  observa¬ 
tion  is  consistent  with  the  loss,  at  a  linear  decay  rate  of  order  machine  epsilon, 
in  the  SMPM  discussed  in  §5.1.  In  contrast,  the  DGM  and  SEM  conserve  mass 
to  the  order  of  machine  epsilon  regardless  of  the  value  of  At.  On  the  other  hand, 
conservation  of  energy  is  improved  by  the  three  methods  once  the  polynomial 
degree  increases  or  the  size  of  At  decreases. 
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a.  h  Normalized  Lg  error  norm 


b.  h  Normalized  L2  error  norm 


Polynomial  degree  (N) 


c.  h  Normalized  L2  error  norm 


Figure  13:  Convergence  plots  for  the  Standing  wave  problem  when  different  At  sizes  are  used. 


a.  Mass  conservation  (m2/s2) 
SMPM 


b.  Mass  conservation  (m2/s2) 
DG 


c.  Mass  conservation  (m2/s2) 
SEM 


Figure  14:  Conservation  of  mass  for  the  Standing  wave  problem  when  different  At  sizes  are 
used. 


5.3.  Effect  of  Filtering 

In  the  interfacial  treatment  of  the  SMPM,  there  is  no  dissipative  term  that 
removes  spurious  high  wavenumber  oscillations  that  develop  in  highly  nonlinear 
simulations.  Thus,  spectral  filtering  is  needed  when  such  simulations  are  run 
for  long  integration  times,  namely  when  sharp  localized  features  emerge  in  the 
simulations  (e.g.,  nonlinear  Riemann  problem)  or  even  when  the  structure  of  the 
solution  is  apparently  smooth  and  free  of  any  localized  features  (e.g.,  nonlinear 
Stommel  problem).  In  contrast,  in  the  case  of  the  DGM  ,  the  dissipation  term 
introduces  a  dissipation  mechanism  that  stabilizes  the  solution  and  renders  it 
oscillation-free;  for  a  very  simple  flow  problem,  this  term  reduces  to  a  simple 
upwinding  scheme.  By  neglecting  it,  the  DGM-generated  solution  also  becomes 
unstable.  Without  resorting  to  recasting  the  nonlinear  terms  in  skew-symmetric 
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Relative  error  of  Energy 


a.  Energy  conservation  (m3/s2) 


b.  Energy  conservation  (m3/s2) 


c.  Energy  conservation  (m3/s2) 


Figure  15:  Conservation  of  energy  for  the  Standing  wave  problem  when  different  At  sizes  are 
used. 


form  [7]  and  in  the  absence  of  an  over-integration-based  de-aliasing  strategy  [53] 
(both  which  are  out  of  the  scope  of  this  paper),  spectral  filtering  is  required  to 
recover  stability.  In  terms  of  mass  and  energy  conservation,  the  performance  of 
the  DGM  appears  to  be  very  similar  when  spectral  filtering  and  no  dissipative 
term  is  used  or  when  only  the  dissipative  term  is  used  (Fig.  12). 

The  performance  of  both  SMPM  and  DGM  is  further  examined  in  problems 
where  significantly  sharp  features  are  present.  The  darn-break  problem  [9]  is 
simulated  with  a  cylindrical  step-function  of  the  free  surface  elevation  as  an  ini¬ 
tial  condition  and  with  (. x,y )  €  [—20,20]  x  [—20,20]  and  t  G  [0,0.1].  The  effect 
of  filtering  (with  a  Boyd-Vandeven  filter  of  p  =  10)  and  the  dissipative  term  on 
the  solution  are  shown  in  Fig.  16. 

In  the  absence  of  a  dissipative  term  in  the  DGM  and  any  spectral  filtering  for 
both  methods  (panels  a  and  d) ,  spurious  oscillations  are  localized  in  the  vicinity 
of  subdomain  interfaces  for  the  SMPM,  whereas,  in  the  DGM,  these  oscillations 
are  more  evenly  distributed  throughout  the  computational  domain.  When  spec¬ 
tral  filtering  is  applied  to  both  methods  (panels  b  and  e),  the  oscillations  are 
strongly  damped  in  the  subdomain  interior  where  the  effect  of  the  filter  is  fo¬ 
cused  [11].  Nevertheless,  some  weaker  oscillations  remain  at  the  subdomain 
interfaces  [11].  If  no  spectral  filtering  is  applied  but  an  additional  dissipative 
term  is  added  to  the  penalty  term  in  the  SMPM  (panel  c),  the  solution  has 
a  near  identical  structure  with  the  one  computed  by  the  DGM  with  the  full 
Rusanov  flux.  For  the  purpose  of  comparison,  Fig.  17  shows  the  filtered  solu¬ 
tion  obtained  from  the  SEM  which  is  contrasted  to  its  filtered  counterparts  (no 
Rusanov  flux  term  present)  computed  from  DGM  and  SMPM  (Figs.  16b  and 
e).  The  results  for  SEM  with  filtering  show  stronger  spurious  oscillations  than 
SMPM  or  DGM  with  dissipation  or  spectral  filtering. 
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a.  SMPM  no  dissipation 


c.  SMPM  with  dissipation 


i  r 


3r  f.  DGM  with  dissipation 


Figure  16:  Cross  section  of  the  Dam-break  problem  for  5x5  subdomains,  and  N  =  20  at 
t  =  0.1.  (a)  SMPM  without  dissipation,  (b)  SMPM  with  filtering  (Filter  order  p  =  10).  (c) 
SMPM  with  dissipative  term,  (d)  DGM  without  dissipative  term  (|A|(qR  —  qL).  (e)  DGM 
with  filtering  (Filter  order  p  =  10),  and  without  dissipation  term,  (f)  DGM  full  Rusanov  flux. 


5-4 ■  Computational  Efficiency  and  implementation 

For  all  test  cases,  the  order  of  magnitude  of  the  CPU  time  per  time  step  has 
been  found  to  be  comparable  for  both  DGM  and  SMPM  and  increases  when 
the  number  of  degrees  of  freedom  increases  due  to  h  or  p  refinement.  Fig¬ 
ure  18a  shows  the  computational  time  for  all  three  methods  considered  in  this 
manuscript  (SMPM,  DGM  and  SEM)  for  different  values  of  N  for  the  Riemann 
problem  with  5x5  subdomains  and  the  same  time  step  value  for  each  method, 
corresponding  to  Courant  Number  =  0.4. 

Figure  18b  shows  the  time  needed  to  advance  a  simulation  to  the  same  final 
time  as  Fig.  18a,  where  the  Courant  Number  is  set  to  the  empirically  computed 
maximum  value  that  enables  a  stable  simulation  for  each  method.  SEM  sim¬ 
ulations  are  found  to  support  double  the  maximum  Courant  Number  value  of 
DGM  and  SMPM  and  are  thus  twice  as  fast.  DGM  and  SMPM  simulations 
were  also  performed  with  a  Courant  Number  value  slightly  above  the  empiri¬ 
cally  obtained  stable  limit  value.  In  this  case,  DGM  was  found  to  destabilize 
faster  than  SMPM. 

Theoretical  justification  for  these  observations  is  gained  by  examining  the 
eigenvalue  spectra  of  the  discretized  1-D  linear  advection  operator  for  each  of 
the  three  discretization  methods  for  a  periodic  domain  with  5  subdomains  and 
N  =  4  (Fig.  19).  In  the  absence  of  the  dissipation  term  in  DGM,  and  as 
expected,  all  three  methods  have  purely  imaginary  eigenvalues.  The  extreme 
eigenvalues  of  DGM  are  roughly  25%  larger  than  their  SMPM  counterparts  and 


30 


Figure  17:  Cross  section  of  the  filtered  Dam-break  problem  for  5x5  subdomains,  and  N  =  20 
at  t  =  0.1  (Filter  order  p  =  10)  for  SEM. 


Figure  18:  CPU  time  for  the  Riemann  problem.  5x5  subdomains  with  different  polynomial 
orders  at  t  =  0.2.  (a)  All  methods  with  Courant  Number  =  0.4  and  (b)  DGM  and  SMPM 
with  Courant  Number  =  0.4  and  SEM  with  Courant  Number  =  0.8. 


double  the  corresponding  SEM  eigenvalues.  Incorporation  of  the  numerical  flux 
term  in  DGM  gives  rise  to  eigenvalues  with  a  negative  real  part  which  equip 
the  numerical  solution  with  the  necessary  numerical  dissipation.  Moreover,  the 
separation  between  the  eigenvalues  with  the  largest  absolute  imaginary  values 
is  reduced  with  respect  to  the  case  without  dissipation  but  is  still  slightly  larger 
than  that  in  SMPM  and  almost  double  that  of  SEM.  Taking  into  account  the 
stability  region  of  the  SSP-RK34  scheme  (which  is  stable  along  the  imaginary 
axis)  for  Courant  numbers  below  this  eigenvalue  separation  can  explain  why 
SEM  can  attain  double  the  Courant  Number  of  DGM  and  SMPM  and  why 
DGM  explodes  a  little  faster  than  SMPM  for  a  marginally  unstable  time  step. 

In  terms  of  implementation,  in  the  context  of  the  SWE,  both  the  SMPM  and 
DGM  can  be  written  as  a  system  of  time-dependent  ordinary  differential  equa¬ 
tions  where  the  vector  of  unknowns  is  the  solution  vector  at  the  grid  points  [54] . 
In  the  matrix-vector  product  that  appears  on  the  right  hand  side  of  this  system 
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Figure  19:  Eigenvalue  distribution  of  the  ID  discrete  linear  advection  operator  (u )  for  all 
three  methods,  with  an  advective  velocity  of  u  =  1.  In  all  cases  x  G  [—1,1],  5  subdomains, 

N  =  4. 

of  equations,  the  associated  matrix  is  simply  a  spectral  differentiation  matrix 
(Eq.  (19)  )  for  the  SMPM  due  to  its  underlying  collocation  method  framework 
with  any  modifications  to  this  matrix  incurred  through  communication  with 
points  on  the  edge  of  the  neighboring  subdomain.  Similar  modifications  on 
account  of  the  numerical  flux  term  enter  the  construction  of  the  corresponding 
right  hand  side  matrix  for  the  DGM,  the  core  of  which  is  built  through  additional 
numerical  integration  and,  therefore,  cost.  This  cost  is,  nevertheless,  offset  over 
the  course  of  a  long  unsteady  simulation.  In  summary,  for  hyperbolic  systems 
of  equations,  the  cost  of  SMPM  and  DGM  are  very  similar.  However,  we  expect 
the  SMPM  to  have  an  advantage  when  elliptic  operators  are  introduced  since  the 
addition  of  a  Laplacian  for  the  SMPM  becomes  simply  a  matter  of  introducing 
a  Laplacian  differentiation  matrix  whereas  in  DGM  either  local  discontinuous 
Galerkin  or  interior  penalty  methods  have  to  be  introduced  [55,  56,  57].  For 
SEM,  the  addition  of  Laplacian  operators  introduces  only  a  slight  cost. 

6.  Conclusions 

The  performance  and  properties  of  two  commonly  used  high-order-accuracy 
element-based  spatial  discretization  methods,  spectral  multidomain  penalty  (SMPM) 
and  discontinuous  Galerkin  (DGM),  are  examined  in  the  framework  of  the  invis- 
cid  shallow  water  equations  (SWE).  Whereas  a  previous  comparison  study  [31] 
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focused  on  one-dimensional  conservation  laws  and  considered  a  modally-based 
Galerkin  formulation  of  SMPM  and  DGM,  this  paper  applies  both  techniques  to 
a  system  of  nonlinear  conservation  equations  and  considers  them  in  the  more  fre¬ 
quently  used  nodal  form,  in  a  collocation  and  Galerkin  formulation,  for  SMPM 
and  DGM,  respectively.  The  two  methods  are  applied  to  a  suite  of  test  cases 
that  are  of  interest  in  oceanic  shallow  water  flow:  three  linear  (standing  wave, 
Kelvin  wave  and  linear  Stommel  problem)  and  three  non-linear  (Rossby  soliton, 
nonlinear  Stommel  problem  and  Riemann  problem).  The  analysis  shows  that 
the  methods  can  be  simplified  to  be  the  same  method  when  specific  choices  of  the 
penalty  terms  (for  the  SMPM)  and  numerical  flux  (for  the  DGM)  and  when  the 
same  collocation  points  are  chosen  for  representing  the  discrete  solution.  The 
numerical  solutions  showed  that  the  methods  are  extremely  similar  not  only  in 
achieving  the  same  rate  of  convergence  but  also  in  their  conservation  of  energy 
measures.  The  key  difference  between  the  SMPM  and  DGM  is  in  their  choice  of 
penalty  terms  that  enforce  weak  boundary  conditions  across  element  interfaces. 
The  SMPM  has  much  flexibility  in  selecting  these  terms  whereas  the  DGM 
method  is  more  rigid  in  its  choices  in  the  sense  that  a  Riemann  solver  must 
be  used;  however,  this  idea  offers  much  flexibility  in  handling  a  large  variety 
of  flows  including  those  requiring  wetting  and  drying  algorithms,  for  example. 
Both  methods  can  be  used  on  fully-unstructured  quadrilateral  element  grids  but 
it  is  not  clear  how  to  extend  the  SMPM  to  unstructured  triangular  elements; 
in  contrast,  the  formulation  of  the  DGM  is  quite  natural  and  can  be  extended 
to  triangles  rather  straightforwardly,  assuming  that  a  good  set  of  interpolation 
and  integration  points  are  known  (see,  e.g.,  [21,  23,  25,  6]).  The  SMPM  proved 
to  be  slightly  more  efficient  than  the  DGM,  in  terms  of  computational  time,  and 
we  expect  this  trend  to  continue  as  Laplacian  operators  (as  required  by  Navier- 
Stokes  or  even  by  more  realistic  shallow  water  ocean  modeling  simulations)  are 
introduced. 

To  observe  further  differences  between  both  methods,  test  cases  with  complex 
geometries,  non-smooth  solutions  or  additional  forcing  terms  have  to  be  exe¬ 
cuted  with  the  methods.  Additionally,  parabolic  and  elliptic  partial  differential 
equations  have  to  be  assessed  in  the  context  of  compressible  and  incompressible 
flows,  where  more  challenging  numerical  difficulties  appear  for  the  implementa¬ 
tion  of  both  methods. 
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